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Abstract. High precision measurements of the Cosmic Microwave Background 
(CMB) anisotropies, as can be expected from the planck satellite, will require high- 

acc'uracy theoretical predictions as well. One possible source of theoretical uncertainty 
is the numerical error in the output of the Boltzmann codes used to calculate angulax 
power spectra. In this work, we carry out an extensive study of the numerical 
accuracy of the public Boltzmann code CAMB, and identify a set of parameters which 
determine the error of its output. We show that at the current default settings, the 
cosmological parameters extracted from data of future experiments like Planck can be 
biased by several tenths of a standard deviation for the six parameters of the standard 
ACDM model, and potentially more seriously for extended models. We perform an 
optimisation procedure that leads the code to achieve sufficient precision while at the 
same time keeping the computation time within reasonable limits. Our conclusion 
is that the contribution of numerical errors to the theoretical uncertainty of model 
predictions is well under control - the main challenges for more accurate calculations 
of CMB spectra will be of an astrophysical nature instead. 
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1. Introduction 

Any meaningful quantitative analysis of experimental data is based on a comparison 
with the predictions of a theoretical model, and the analysis of Cosmic Microwave 
Background (CMB) anisotropies is no exception to this rule. The amount of information 
to be gained from observations is limited not only by experimental uncertainties (such as 
detector noise, foregrounds, etc.), but also by the abihty to accurately predict observable 
quantities from a given theory. There are various sources of theoretical uncertainties. 
Some, such as cosmic variance, are endemic to the problem, and unavoidable. Others 
are based on insufficient theoretical understanding of the complex processes involved [1] 
(examples include the physics of recombination [2], reionisation [3], and contributions 
due to the Sunyaev-Zel'dovich effects). Additionally, the analysis may be compromised 
by inaccuracies in the numerical programmes used to calculate the anisotropy angular 
power spectra [4]. 

With increasingly sophisticated experiments, and the contribution from experimen- 
tal errors becoming less and less important, the relative significance of theoretical errors 
increases. In fact, for the PLANCK satellite [5], the signal uncertainty of the temperature 
anisotropies will be dominated by cosmic variance instead of noise over a wide range of 
scales up to multipolcs of £ ~ 2000. 

Ignoring any type of uncertainty can lead to biased estimates of parameter values, 
and, in the worst case, a wrong physical interpretation of the data. It is therefore 
imperative to either make sure that the errors are small enough to be negligible, or, if 
that should not be the case, to devise an appropriate strategy to deal with the problem. 

In the present work, we shall consider the numerical accuracy of the Boltzmann 
codes employed to calculate the angular anisotropy spectra for given input values of 
cosmological parameters. The first public Boltzmann code was released more than a 
decade ago [6], and to date, there are several other such programmes freely available 
for download [7-9]. The output of all these codes is necessarily inaccurate to some 
extent, due to the use of semi-analytical approximations as well as artifacts of the 
numerical implementation, such as finite integration steps or the need to interpolate. 
These effects can be parameterised by a set of accuracy parameters, whose settings 
determine the accuracy of the output, but also the computation time. Here, we will 
focus our analysis on the CAMB code by Lewis, Challinor and Lasenby [8] in order to 
avoid possible systematic effects caused by differences between codes - a comparison of 
different (more or less) independent Boltzmann codes was performed by Seljak et al. [4], 
who found an excellent qualitative agreement. | We extend their line of reasoning and 
present a detailed analysis of the potential effects of numerical inaccuracies on parameter 
estimates from PLANCK data. Additionally, we optimise the accuracy settings of CAMB 
to find an ideal balance between precision and execution time. 

We proceed in Section 2 by defining an appropriate measure of accuracy, identifying 

I We verified that after updating various parts of CMBf ast (values of pliysical constants, recombination 
code, etc.), the output of CAMB and CMBf ast agree sufficiently well. 
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Table 1. Free cosmological parameters of the model, fiducial values used to generate 
the mock data, and prior ranges adopted in the analysis of Section 4. 



Parameter Fiducial Value Prior Range 



Dark matter density 




0.10976 


0.01 - 


^ 0.99 


Baryon density 




0.02303 


0.005 


0.1 


Hubble parameter 


h 


0.7 


0.4^ 


1 


Optical depth to reionisation 


T 


0.09 


0.01 - 


^ 0.5 


Scalar spectral index 


ns 


0.96 


0.5 ^ 


1.5 


Amplitude of scalar spectrum i 


i k = 0.05 Mpc-i In [lO^°As] 


3.135 


3^4 





the relevant parameters which affect the accuracy of the output spectra and constructing 
a set of CMB reference spectra. In Section 3, we will describe our optimisation procedure 
and present a recommended set of accuracy parameters for CAMB, followed by an analysis 
of the expected potential bias on the cosmological parameters of the vanilla model caused 
by numerical inaccuracies in Section 4. The impatient reader may prefer to skip directly 
to our conclusions in Section 5. 

2. Reference spectra 

2.1. The fiducial model 

In principle, the impact of individual accuracy parameters and the number of relevant 
parameters will depend on the underlying cosmological model assumed, and on the 
values of the cosmological parameters. In this analysis, we will stick to the 6-parameter 
ACDM- "vanilla" -model and we limit ourselves to a point in the space of cosmological 
parameters that lies close to the best fit to the WMAP 5-year data [10], see Table 1. 
We neglect the effects of weak gravitational lensing on the CMB spectra [11] (which 
will of course have to be taken into account in an analysis of real PLANCK data), and 
ignore signatures of non-minimal models, like massive neutrinos, tensor modes, spatial 
curvature, etc., and defer their treatment to future work. 

2.2. Measuring accuracy 

In order to quantify the accuracy of the Boltzmann code output C^"*, we require two 
things: a reference point Q''^ to compare with, and a measure of accuracy. 

The reference spectra C}^^ would ideally be the exact prediction of the theory. In 
practice however, we have to make do with getting close enough to these ideal spectra. 
We will return to this issue and describe the construction of C}*^^ in 2.4. 

To measure accuracy, one might be tempted at first glance to look at the relative 
difference of the spectra (C°"* — for each £. However, this approach does 

not properly take into account the fact that the accuracy requirements are dependent 
on i (due to cosmic variance and experimental errors). Additionally, the spectra are 
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merely an intermediate step in the inference process. In the end, we are interested in 
possible biases on cosmological parameters, not the accuracy of the spectra. A more 
meaningful measure of deviation from the reference spectra is the effective x^, which 
can be obtained by using the reference spectra to generate a fiducial data set, taking 
into account the experimental errors of an experiment (or the projected errors in case 
of future experiments) , and "fitting" Q to these data. 

More formally, the effective is related to the likelihood C and defined by 

Tr(C7^Q)+ln 



X^ = -21n£ = ^(2^+l) 



- 2 



(2.1) 



Here, is the theoretical covariance matrix, and its entries are taken to be the sum of 
the signal and noise power spectra: 

C,^{cr+Nr'], (2.2) 
where the index X runs over temperature (T) and polarisation (£■), and for a fiducial 
data set, we can take the data covariance matrix to be equal to , i.e., the 

do 

theoretical covariance matrix evaluated for the fiducial values of the cosmological 
parameters. 

We take the noise to be isotropic and Gaussian; the noise power spectrum is related 
to the experimental parameters of Table 2 through 



XX' 



SXX' ^beam^X exp 



£{£+![ 



n2 

beam 

8 In 2 



(2.3) 



For more details see Refs. [12, 13]. It should be noted that the normalisation of Eq. (2.1) 
is chosen such that the total is zero when the output spectra exactly match the 
reference spectra used to construct the fiducial data. 

We generate the fiducial data set of TT-, EE-, and T£'-spectra up to £ = 3000 
using the code of Perotto et al. [12] . For simphcity we ignore the effects of incomplete 
sky coverage due to masking the galaxy and point sources, as well as anisotropic noise 
[14]. To evaluate the accuracy of CAMB in view of the expected data from PLANCK, we 
assume 14 months of integrated observations in the 70 GHz channel of LFI and the 100 
and 143 GHz channels of the HFI instrument; their specifications are taken from the 
PLANCK blue book [5] and listed in Table 2. 



2.2.1. Interpretation of the measure As can be seen in Eq. 2.1, is directly 
related to the likelihood, which, along with a choice of prior probability densities on 
all cosmological parameters, leads to the posterior probability density, from which 
constraints on parameters are eventually derived. Assuming fiat priors on the 
parameters and a multivariate Gaussian likelihood function, for a given numerical error 
X^, the bias on any cosmological parameter cannot exceed VX^ standard deviations in 
the worst case (i.e., when the error in the angular power spectra can be exactly offset by 
changing one of the cosmological parameters). On the other hand, if the error had no 
degeneracy with any cosmological parameter, the effect of a non-zero x^ would be just a 
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Table 2. List of technical specifications for the 70 GHz channel of the LFI and the 100 

and 143 GHz channels of the HFI instrument: 6'beam denotes the beam width, At,p 
are the sensitivities per pixel and u is the centre frequency of the channels. 



i//GHz 


^bcam 




Ap/^K 


70 


14.0' 


12.8 


18.3 


100 


9.5' 


6.8 


10.9 


143 


7.1' 


6.0 


11.4 



constant offset to the likelihood function, which would not have any effect on parameter 
inference. In realistic cases, the expected bias would lie somewhere in between. For the 
parameters of the vanilla model, we will provide an estimate of the bias in Section 4. 

2.3. The accuracy parameters 

The numerical accuracy of a Boltzmann code's output depends on many factors, 
ranging from the use of analytical approximations to the sampling various intermediate 
quantities in the calculation. These sources of numerical error can be quantified in terms 
of accuracy parameters, e.g., the number of samples used for interpolating a particular 
quantity. An increase in accuracy will generally be accompanied by a longer computation 
time and possibly higher requirements on the available computer memory. 

We use the June 2008 version of CAMB§ as a starting point of our analysis. The 
unmodified version of CAMB comes with a set of three continuously adjustable accuracy 
parameters: [| 

• l_sample_boost: determines for which values of i. the Ci are actually calculated 

(the rest are interpolated). 

• l_accuracy_boost: determines the multipole at which the Boltzmann hierarchy for 
photons, neutrinos, etc., is cut off. 

• accuracy.boost: affects the setting of various time steps, samplings, etc. 

The latter two parameters affect several settings at once, so in the interest of optimising 
the performance of CAMB, we split them up into their constituents and treat them 
separately. Apart from the settings governed by these three parameters, we identified a 
few other quantities which can affect the accuracy of the results, and should be taken 
into account when optimising the code. Altogether, we consider a set of 19 accuracy 
parameters in our analysis, listed in Table 3. 

All parameters are defined in such a way that setting them a value of 1 reproduces 
the results of the unmodified version of CAMB, and larger values correspond to 
better accuracy. The constituents of the old 1 accuracy boost and accuracy boost 

§ http: //www. camb. info 

II There are also a few logical switches that arc relevant here; in our analysis wc kept them fixed to 
accurate_polarization=T, accurate_reionization=T and do_late_rad_truncation=T. 
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Parameter 


Corresponding old accuracy parameter 


Description and comments 


new_l_sample_boost 


l_saiiiple_boost 


^-sampling of C^s 


Imaxg 
Imaxnr 


l_accuracy_boost 


Boltzmann hierarchy cutoff for photons 
Boltzmann hierarchy cutoff for massless neutrinos 


int_tol 

i _'t i m e s "t e p 
rec_timestep 
rec_tiinestep2 

rad_trunc 
deC-Start 

-i- ±± V G W-iil k/ -1- w 

source_dk 
source_kiniii 
int-xlmaxl 
tc_largek 


accuracy_boost 


tolerance parameter for integration routines 

time step during reionisation 

time step during recombination 

time step between recombination and reionisation 

truncation of photon hierarchy during matter domination 

starting time of decoupling 

samnles for inteffration over source function 

fc-sampling of source function 

minimum value of k to calculate source function for 
starting time for source function integration 
switch off tight coupling later for large k 


tc_epO 




tight coupling switch 


bess-sampling 

bess-xlimmin 

bess-xlimfrac 




x-sampling of spherical Bessel functions je{x) 
approximate je{x) ~ for small a;, if ^ > xlimmin 
approximate je{x) ~ for large if a; < (1 — xlimf rac) • i 


ketamajc 




maximum value of kr] 



parameters are taken to multiply the old l_accuracy_boost and accuracy_boost 
parameters (e.g., setting Imaxg = Imaxnr = 2 produces the same effect as setting 
l_accuracy_boost = 2). We modified the routine that determines for which values of 
I the C(, are calculated: our parameter new_l_sample_boost is defined to be the square 
root of the old l_sample_boost; for new_l_sainple_boost > 5, all Ce, are calculated and 
there will be no interpolation of the final spectrum. 

2.4.- Constructing a reference data set 

To quantify absolute accuracy we require a reference data set, as discussed above. Its 
construction is naturally tied to choosing the accuracy parameters in such a way that 
increasing them further would not have any appreciable effect. However, by arbitrarily 
increasing all parameters to "large" values at the same time, one would run into the 
limits of the hardware, particularly the available memory. For the purpose of finding 
suitable values for generating the reference spectra, we therefore analyse the parameters 
one by one, keeping all other parameters fixed. For each parameter, we generate a 
fiducial reference data set with that parameter set to a high value, all other parameters 
kept at a value of 2. Varying this parameter and calculating the reveals its impact 
on overall accuracy, allows us to find a suitable setting for the reference spectra and lets 
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us estimate the remaining error. 

We show the results of the single parameter scans in Figure 1-2. From Figure 1 
we can see that certain parameters (e.g., Imaxg, Imaxnr) are very well-behaved and 
reach a ~ 10~® plateau well below the fiducial value. For these parameters it is 
reasonable to assume that picking even higher values would not have any appreciable 
effect on accuracy. A number of other parameters do not fully converge, and exhibit 
a step-like behaviour before reaching the fiducial value (e.g., int .sample, source.dk). 
It is likely that increasing their value beyond our fiducial maximum value would have 
an effect on the spectra. However, the graphs in Figure 1 nonetheless provide an order 
of magnitude estimate of the remaining error and one can use them as a guide to 
finding reasonable settings for the construction of the reference spectrum. Finally, 
the parameter source_kmin does not seem to converge at all and exhibits unstable 
behaviour, though its effect on overall accuracy is negligible. 

The parameter settings we chose for the reference spectrum are given in Table 4. 
The dominant contribution to any residual error of the reference spectrum will come 
from the parameter displaying the worst convergence - ketamax. Unfortunately, this 
parameter also has a strong impact on the computation time T (see Figure 2), and 
memory requirements, precluding us from choosing a higher setting. We estimate the 
reference spectrum to be accurate to Ax^ of order 10~^. 

3. Optimising performance 

Having constructed the reference spectra, we can now proceed to searching settings of 
the accuracy parameters which yield a result that hes as close as possible to the reference 
spectra, within a reasonable time of computation. To this end, we use a modified version 
of the Markov chain Monte Carlo code CosmoMC [15]. Major modifications include: 

• we use the fiducial data set constructed from the reference spectra; 

• we vary the accuracy parameters instead of the cosmological parameters (which are 
kept fixed at their fiducial values) , taking top hat priors on all accuracy parameters 
(with lower limits at a value of 0.5 and upper limits large enough to not infiuence 

the results); 

• instead of sampling from the usual posterior V (which is proportional to 
C ~ exp[— x^/2]) itself, we sample a function F{V), defined by 



This function was chosen such that areas of parameter space leading to too long 
computation times are avoided, and that areas of parameter space giving <^ 1 
are better sampled. 

We generated ~ 20000 sample settings in this way; a scatter plot of the samples in 
the (x^,T)-plane is presented in Figure 3, illustrating the strong correlation between 



FiV) = 



if T > 60 s 

V if P > 1 

_ 5 In [P] + 1 if V <l 



(3.1) 
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Table 4. Accuracy settings for the reference spectrum (plus three recommended 
example settings, see Section 3). 



Parameter 


Reference setting 


Setting 1 


Setting 2 


Setting 3 


new_l_sample_boost 


6 


1.77 


1.60 


2.00 


Imaxg 





T on 
7.39 


z.Oz 


one 

z.Oo 


Imaxnr 





1.78 


Z.ll 


O.OD 


int_tol 


6 


2.07 


1.53 


5.86 


ri_timestep 


2 


0.99 


0.71 


2.87 


rec_timestep 


3 


0.50 


0.50 


0.75 


rec_timestep2 


10 


1.85 


1.14 


2.44 


rad_trunc 


4 


1.81 


1.76 


2.58 


dec_start 


100 


2.35 


35.76 


5.05 


int .sample 


6 


4.21 


4.88 


3.65 


source.dk 


4 


3.10 


2.84 


2.72 


source_kmin 


1 


5.86 


2.50 


5.23 


int_xlmaxl 


4 


1.20 


1.00 


2.16 


tc_largek 


5 


2.44 


2.33 


1.90 


tc_epO 


5 


4.25 


3.30 


6.32 


bess_sampling 


3 


3.02 


2.53 


4.00 


bess_xliminin 


2 


1.14 


0.74 


3.61 


bessjxlimf rac 


2 


4.19 


0.90 


1.07 


ketamax 


3 


1.32 


0.96 


0.56 



accuracy and computation time. Which settings to use is a somewhat subjective 
decision and should be taken with one's available computing power in mind. We list two 
sample settings in Table 4, taken from the samples of our Monte Carlo run: "Setting 1" 
corresponds to the best accuracy under the constraint that the computing time be less 
than 60 seconds, "Setting 2" gives the best accuracy for T < 30 s, and "Setting 3" the 
best accuracy for T < 17.6 s (which is the computation time for the unmodified version 
of CAMB run at accuracy_level = 2). The performance of these settings is contrasted 
to game's default (accuracy.level = 1), and a high-accuracy setting of an unmodified 
version of CAMB (accuracy_level = 2) in Table 5. 

Note that at the default settings, with a difference of Ax^ to the reference spectra, 
parameter estimates could be biased by more than two standard deviations, in the worst 
case. For both Settings 1 and 2, on the other hand, the maximum possible bias would be 
less than 0.1 standard deviations (~ 0.13 standard deviations for Setting 3), assuming 
Gaussian posterior distributions. 

The results of the MCMC search confirm the tendencies observed in the single 
parameter scans of Figs. 1 and 2 regarding the impact of individual parameters on 
accuracy and speed; we find no evidence for significant cross-correlations between 
accuracy parameters. 



Boltzmann code optimisation 



9 



Table 5. Deviation from reference spectrum, given in terms of x^j and computation 
time T, for the three example settings of Table 4, CAMB's default settings 
(accuracy_level = 1) and accuracy .level = 2. 



Setting 










1 






2.7- 10-3 


59 


2 






7.0 • 10-3 


27 


3 






1.6-10-2 


17.4 


accuracy. 


-level = 


1 


5.8 


2.7 


accuracy- 


-level = 


2 


0.5 


17.6 



4. Estimating the bias 

As mentioned above, the for a given accuracy settings can be used to estimate the 
worst case bias on cosmological parameters. Nevertheless, we are also interested in how 
large a bias effect we can expect in the example case of the vanilla model. To this end, 
we performed an actual parameter estimation exercise using CosmoMC, for three cases: 

• Using the same CAMB accuracy settings for generating the fiducial data as for the 
inference process - mimicking an analysis free of numerical errors. 

• Using the reference data set and running CAMB at accuracy settings "2" from Table 4. 

• Using the reference data set and running CAMB at its default settings. 

We generate 16 Markov chains, making sure the Gelman-Rubin convergence parameter 
R — 1 [16] is smaller than IQ-^ for all parameters considered. The results are plotted 
in Figure 4: as expected, there is no discernible bias between the results obtained with 
optimised settings and an error-free analysis. The expected bias for the default settings 
is rather mild, not exceeding a few tenths of the posterior standard deviations for the 
respective parameters, with ns (40%) and In [10^° As] (43%) being the most affected. 

5. Conclusions 

To meet the challenge posed by ultra-precise future CMB experiments, the output of 
the Boltzmann codes used to calculate theoretical predictions for the CMB anisotropy 
spectra will have to meet stringent requirements in terms of numerical accuracy - a 
goal that comes at the cost of accordingly higher demand in computational resources. 
This tendency is partially mitigated by the increase in available computing power. 
Additionally, the use of efficient interpolation algorithms [17-20] can significantly 
accelerate the process of parameter estimation from future data sets. Nonetheless, 
such methods still rely on the input of a full Boltzmann code, whose inherent numerical 
errors will be propagated to the interpolation codes. 

We have performed a detailed analysis of the numerical accuracy of CAMB, provided 
an estimate of the residual numerical error of the output, and evaluated the possible 
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impact on parameter estimates from PLANCK data. For the present default accuracy 
settings of CAMB and simulated PLANCK data, numerical errors can lead to a slight bias 
on estimates of the six free parameters of the ACDM model; though we cannot exclude 
the possibility of a more serious bias on parameters of extended models. However, by 
tweaking the settings of various internal parameters of CAMB it can be made sure that 
the bias on any parameter will not exceed 0.1 standard deviations. An example for 
suggested settings of the accuracy parameters is given in Table 4. 

The results of this paper lead to an efficient and more accurate calculation of 
CMB angular power spectra, and should bring CAMB to a standard that will allow us 
to make the most out of upcoming PLANCK data. We conclude that the contribution 
of numerical errors to the theoretical uncertainty of model predictions is well under 
control - the main challenges for more accurate calculations of CMB spectra will be of 
an astrophysical nature instead. 
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Figure 1. These diagrams illustrate the dependence of accuracy on the settings of 
individual parameters. We plot ^ ^ function of parameter value. The fiducial 
reference data sets were generated with all parameters set to a value of 2, except for 
the one scanned, which is set to an extremely high value (100 for dec_start, 5 for 
tc_largek, and 10 for all other parameters). When setting the accuracy parameters 
to their fiducial values, we obtain a of order 10"* instead of 0. This effect stems 
from a numerical rounding error when outputting the fiducial data set, and is small 
enough not to be of relevance to any of our conclusions. 
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X 



X 
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Figure 2. These diagrams show the computation time T (in seconds), varying the 
settings of all parameters individually while keeping all others fixed at a value of 2. 
The calculations were performed on a single core of a 2.40 GHz Intel T7700 CPU. The 
spikes for new_l_sample_boost, as well as the "bumps" at low settings of various other 
parameters are due to the fact that the Bessel functions were (re-)calculated at these 
points, adding a few seconds to the total time. 
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Figure 4. Marginalised posterior probability densities of the cosmological parameters 
of the vanilla model. Thick black lines correspond to the exact results (using the 
same accuracy settings for the MCMC and for the fiducial data set), green lines are 
results obtained with the "Settings 2" column of Table 4 and the reference data set, 
and the red dotted lines represent the posteriors from the default accuracy settings of 
CosmoMC/CAMB and the reference data set. 





